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Abstract 

We study the numerical performance of a continuous data assimi¬ 
lation (downscaling) algorithm, based on ideas from feedback control 
theory, in the context of the two-dimensional incompressible Navier- 
Stokes equations. Our model problem is to recover an unknown ref¬ 
erence solution, asymptotically in time, by using continuous-in-time 
coarse-mesh nodal-point observational measurements of the velocity 
field of this reference solution (subsampling), as might be measured 
by an array of weather vane anemometers. Our calculations show 
that the required nodal observation density is remarkably less that 
what is suggested by the analytical study; and is in fact comparable 
to the number of numerically determining Fourier modes, which was 
reported in an earlier computational study by the authors. Thus, this 
method is computationally efficient and performs far better than the 
analytical estimates suggest. 
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1 Introduction 


The goal of data assimilation is to provide a more accurate representation 
of the current state of a dynamical system by combining observational data 
with model dynamics. This allows the influences of new data to be incor¬ 
porated into a numeric computation over time. Data assimilation is widely 
used in the climate sciences, including weather forecasting, environmental 
forecasting and hydrological forecasting. Additional information and histor¬ 
ical background may be found in Kalnay [13] and references therein. 

In 1969 Charney, Halem and Jastrow [5] proposed a method of continuous 
data assimilation in which observational measurements are directly inserted 
into the mathematical model as it is being integrated in time. To fix ideas, 
let us suppose that the evolution of u is governed by the dynamical system 

du 

— = A(m), M(to) = (1) 

and the observations of u are given by the time series p{t) = Pu{t) for 
t G [to,t^], where P is an orthogonal projection onto the low modes. In this 
context, the method proposed in [S] for approximating u from the observa¬ 
tional data is to solve for the high modes 


= {I - + (l{to) = qo ( 2 ) 

where qo is an arbitrarily chosen initial condition and q + p represents the 
resulting approximation of u. Note that if qo = {I — P)uq then p + q = u ioi 
all time; however, data assimilation is applied when uo is not known. 

Algorithm (|2]) was studied by Olson and Titi in mi and [18] for the 
two-dimensional incompressible Navier-Stokes equations 


du 

'm 


uAu -I- (u ■ V)m -I- Vp = / 
V • M = 0 


(3) 


on the domain D = [0,L]^, equipped with periodic boundary conditions and 
zero spatial average with initial condition u{x,to) = uo{x) for x G [0,L]^. 
Observational measurements were represented hj P = Ph, where Ph is the 
orthogonal projection onto the Fourier modes exp{27rik ■ x/L) with wave 
numbers fc G \ {0} such that 0 < |fc| < L/h. Here i/ > 0 is the kinematic 
viscosity, p{x,t) is the pressure and /(x) is a time-independent body force 
with zero spatial average acting on the fluid. For simplicity, it was assumed, 
as we shall here, that V ■ / = 0. 
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The two-dimensional incompressible Navier-Stokes equations are amen¬ 
able to mathematical analysis while at the same time they posses non-linear 
dynamics similar to the partial differential equations that govern realistic 
physical phenomenon. Using the functional notation of Constantin and Foias 
[6], see also Temam [20] or Robinson [19], write ([3]) in the form ([T]) by setting 

J-{u) = —uAu — B{u,u) + f (4) 

where A and B are the continuous extensions of the operators given by 


A = —P(jAu and B{u,v) = P^iu-Vv) 


when M, V are smooth divergence-free L-periodic functions, and P^ is the 
Leray-Helmholtz projector. We recall that 


E I 

fcGZ2\{0} ^ 


Uk - \ exp{2Trik ■ x/L) 


and also that A-. —)■ V~^ and B: x —)■ V~^ where U" is the closure 
of V, the space of zero-average R^-valued divergence-free L-periodic trigono¬ 
metric polynomials, with respect to the norm 


uWvc^ = \k\^°‘\uk\^. 

fcGZ2\{0} 


For notational convenience we shall write V = throughout the remainder 
of this paper. 

Consider the data assimilation method given by ([2]). Using the theory of 
determining modes it was shown as Theorem 1.5 in im that if h satisfies 

> ciG (5) 

then \\u{t) — pit) — q{t)\\y —)■ 0, exponentially fast as t —)■ cxd. Here Ci 
is a universal constant and G = (L/27rz/)^||/||2,2 is the Grashof number. 
Computations in [TS] considered a hxed forcing function / = /121 scaled 
to obtain different values for G. In that work the subscript 121 was used 
to indicate that the force was supported on an annulus around = 121 in 
Fourier space. More details on / are provided by equation (fTT|) and Figure [T] 
below. For Grashof numbers between 500 000 and 60 000 000 it was shown 
that the projection P^ onto the lowest 80 Fourier modes was necessary and 
sufficient to ensure numerically that ||M(t) — pit) — qit)\\Y —)■ 0, as t —)■ cxd. 
Since the rank of Ph scales as L^/h^, this is signihcantly less than the millions 
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suggested by the analytical bound (jS]). Thus, the data assimilation algorithm 
given by (|2]) performs far better than the analysis suggests. Note, however, 
that this algorithm is not suitable when the observations are given by nodal 
measurements of the velocity held. 

An approach to data assimilation, for dynamics governed by equations 
(151) or equivalently (|1]), which is applicable to nodal observations, was intro¬ 
duced and analyzed by Azouani, Olson and Titi in [T]. Let be a general 
interpolant observable satisfying the approximation identity inequality 

||m - Ih{u)\\l2 < 7i^^II«IIhi + 'l2h^\\u\\H2- ( 6 ) 

Given Ih{u{t)) for t e [to)G], solve 
dv 

— = T{v) + nP„{Ih{u) - v{to)=Vo, ( 7 ) 

where Vq is an arbitrary initial condition. The constant /i is a relaxation 
(nudging) parameter which controls the strength of the feedback control 
(nudging term). In particular, the nudging term pushes the large spatial 
scales of the approximating solution v toward those of the reference solution 
u while the viscosity stabilizes and dissipates the hne spatial scales and any 
spillover into the hne scales caused by the nudging term. It follows from 
Theorem 2 equation (39) of pQ that if h and /i satisfy 

V 

then \\u(t) — v(t)\\v —)• 0 as f —)■ oo. Here C 2 is a universal constant and cq is 
a constant depending only on 71 and 72 of ([ 6 ]). 

The number of nodal measurements needed to uniquely determine a solu¬ 
tion to the two-dimensional Navier-Stokes equations, as f —)■ 00 , was found 
by Foias and Temam in [9] and further rehned in Jones and Titi m Up 
to a logarithmic correction, the analytic bounds given by ([ 8 ]) are the same 
as those given in na. In this paper we check the numerical performance of 
the data assimilation algorithm ([7]) using nodal measurements given by Ih 
for the same body forcing / considered in [15] scaled so that G = 2 500 000. 

Let Qi be disjoint squares that cover [0,L]^ with centers Xi and sides of 
length h = L/K, where K'^ = N. An interpolant operator based on the nodal 
measurements u{xi, f), for i = 1 , 2 ,..., A^, and t G [to, G], which satishes ([ 6 ]) 
is 

Ih{u){x,t) =Ih{u){x,t) - — f Ih{u){x,t)dx, (9) 

J[0,L]2 

where 

N 

Xhiu){x,t) = '^u{Xi,t)xQi{x). ( 10 ) 

i=l 
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We also consider the smoothed interpolation 


Ih{u){x,t) =Xh{u){x,t) - ^ f 

J[o,i 

where 

N 

Xh{u){x,t) = '^u{xi,t){p^ 

i=l 

and Pe{x) = e~‘^p{x/e) with 

P«) = |^““p(r^+r^) 

lo 

and 

To make the smoothing scale compatible with the resolntion parameter 
we take e = ph for some hxed p > 0. An analysis of a smoothed interpolant 
similar to Ih appears in Appendix A of [1] for r; = 0.1 and shows that Ih 
satishes ([6]). Note also that J/i —)■ J/i as r/ —)■ 0. When p is between 0 and 1 
the convolntion rednces the high-freqnencies that wonld otherwise be present 
in the Fonrier series representation of the characteristic fnnction xq.. Valnes 
of p greater than 1 blnr nearby nodal measnrements together. This farther 
downscaling conld be nsefnl in the presence of noisy measnrements, see for 
example [3]; however, the measnrements stndied here will be error free. In 
this work we vary p between 0.1 and 2.0 and hnd that p = 0.7 leads to near 
optimal performance for onr data assimilation experiments. 

In particnlar, onr resnlts show that \\u(t) — v(t)\\v ^ 0, as f ^ oo, when 
the resolntion K of the observational measnrements satishes K > 8 and p 
and p are appropriately chosen. Moreover, if 77 > 9 and p = 0.7, there is 
a wide range of valnes for p snch that the algorithm works well. Since 64 
and 81 nodes are comparable in resolntion to 80 Fonrier modes, the nnmer- 
ical efficiency of algorithm ([7]), nsing nodal measnrements, is comparable to 
algorithm (|2]) nsing Fonrier modes 

Rigorons mathematical analysis of the method of data assimilation stnd¬ 
ied compntationally in this paper has recently been generalized to Benard 
convection by Farhat, Jolly and Titi [7] where it was shown that only ob¬ 
servational measnrements of the velocity held is snfficient to recover the fnll 
reference solntion, i.e., the velocity held and the temperatnre. Inspired by [7j 
Farhat, Lnnasin and Titi [8] have recently improved the algorithm stndied 
here, i.e. the one introdnced in |T], by showing that it is snfficient to nse 


Xh{u){x,t)dx, (11) 


^XQi){x), 


( 12 ) 


for 161, 161 < 1 
otherwise. 


(13) 


di2dii. 
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observational measurements of only one component of the velocity field to 
recover the full reference solution. Further implementation of this algorithm, 
for the subcritical surface quasi-geostrophic equation, has recently been es¬ 
tablished by Jolly, Martinez and Titi HU. The algorithm studied here is also 
closely related the 3DVAR data assimilation method developed by Blomker, 
Law, Stuart and Zygalakis [1] for the Navier-Stokes equations and by Law, 
Shukla and Stuart [lU for Lorenz equations. 

This paper is organized as follows. Section 2 describes the physical pa¬ 
rameters, forcing and initial conditions used to generate the reference solution 
to the two-dimensional incompressible Navier-Stokes equations that we will 
be observing through nodal measurements of the velocity held. Section 3 
reports our computational results, and section 4 gives details of our compu¬ 
tational methods. The last section concludes that data assimilation of nodal 
measurements, by means of equation o. as studied in this paper works 
computationally just as efficiently as equation ([2]) with Fourier modes. 

2 The Reference Solution 

To focus on how the smoothing and resolution of the observational mea¬ 
surements affect algorithm ((U), as well as how to optimize the value of the 
relaxation (nudging) parameter /r, we hx the viscosity and the size of the 
periodic box so that 


V = 0.0001 and L = 27r 

for the remainder of this paper. We further perform all our simulations us¬ 
ing the same reference solution u{t) to the two-dimensional incompressible 
Navier-Stokes equations. As shown in czi the spatial structure of the func¬ 
tion / used to force the reference solution can have a signihcant effect on 
data assimilation. Therefore, to allow comparison with previous results, we 
use the exact same forcing function dehned in El. and further studied in 
[IB] , scaled such that G = 2 500 000 in our present computations. 

This function / is supported on the annulus in Fourier space with wave 
numbers k such that 110 < < 132. In particular, 

f{x)= ^ /fcexp(iA;-a:) (14) 

110<fc2<i32 

with fk = f-k and k ■ fk = 0, where the values of fk are given by Table 1 
in [18] scaled to obtain the desired Grashof number. Note that this forcing 
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is at length scales of about 1/11-th the size of the periodic box. This fact is 
further reflected in the level curves of 

curl / = curl(/i, /a) = ^ ^ 

ox oy 

depicted in Figure [Don the left. 

Figure 1: Left are contours of curl/; right are vorticity contours of Uq. 



The initial condition uq used for our data assimilation experiments was 
obtained by solving (|3]) with zero initial condition at time t = —25000 until 
time t = 0. In terms of eddy turnover times, this ensures that more than 
500 eddy turnovers have occurred before reaching t = 0. Integrating for this 
length of time ensures the initial condition uq lies close to the global attractor 
and therefore reflects the energetics of the forcing /. In particular, the way 
in which we initialized the solution at time t = —25000 is unimportant. 

The vorticity contours uq = curl uq of the initial condition uq are depicted 
in Figure [T] on the right. While the forcing / contains no Fourier modes 
with wave numbers k such that \k\ < 10, the initial condition uq clearly 
possesses two large eddies the size of the box. These large box-hlling eddies 
apparently result from the inverse cascade of energy in the two-dimensional 
Navier-Stokes equations. This can be seen more clearly by examining the 
energy spectrum of the reference solution. 

Given a solution u{t) to the two-dimensional Navier-Stokes equations for 
t G [0,T] where T = 25 000 dehne the average energy spectrum as 

E{r) = — / \uk(t)\^dt where u(t) = ^ 

^ keJr fcGZ2\{0} 

and jLr = { fc G ; r — 0.5 < |A:| < r -|- 0.5 }. For the reference solution 
described above with initial condition uo and time to = 0, the average energy 
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spectrum appears in Figure [2J While we do not see the Kraichnan scaling 
of k~^ in the inertial range, we do see the Kolmogorov scaling in the 

inverse cascade. Such spectra have been observed in other numerical experi¬ 
ments, see, for example, Xiao, Wan, Chen and Eyink 1211 . In particular, the 
inverse cascade appears responsible for the box-hlling eddies observed in the 
initial condition uq which persist for all times t > 0. 

Figure 2: Time averaged energy spectrum of the reference solution. 



We compute the eddy turnover time for the reference solution as 
r = ^ ^ E(r)^ ^30.8 

r=l r=l 

and conclude that our averages have been computed over T/t ^ 812 eddy 
turnovers. The spectrum of r/, which also has units of energy, has been 
plotted in Figure [2] as three circles to illustrate where the forcing lies in 
relation to the energy spectrum. Note that the forcing exactly divides the 
energy spectrum between the part which scales as k~^^^ and the part which 
scales as k~^. 

Having, to some extent, described the reference solution that will be used 
in our numerical experiments, we now turn to our main point of study, the 
data assimilation of nodal measurements of the velocity held. 






3 Nodal Observations of Velocity 

We consider nodal observations u{xi, t), for z = 1,. .., iV, of the reference so- 
Intion M, that was compnted according to the incompressible two-dimensional 
Navier-Stokes eqnations ([3]) and initialized with uq as described in FignreH], 
and interpolate these measnrements nsing the operator Ih dehned by (jH]). 
The resnlting eqnations for the approximating solntion v may be written as 

ulJ 

— + vAv B{v, v) = f - nP^{Ih{u) - Ih{v)) (15) 

where v is initialized as vq = 0. Note that only the observations Ih(u) of 
the reference solntion u enter into the eqnations for compnting for v. Also 
note that ||m( 0) — n(0)||y = UmoIIv ~ 1.946. Onr goal now is to choose the 
resolntion parameter h and the relaxation (nndging) parameter fi in snch a 
way that \\u(t) — v(t)\\v 0, nnmerically, as t —>■ cx). 

As discnssed in [1], if ft is too small, the feedback control (nndging term) 
will be too weak to ensnre the approximating solntion converges to the ref¬ 
erence solntion. If /i is too large, then spill over into the fine scales becomes 
signihcant and again prevents recovery of the reference solntion. FignreOil- 
Instrates each of these possibilities for h = L/K, where K = 9, nsing different 
valnes of fi. 

Fignre 3: The error ||M(t) — 'i^(t)||y versns t for h = 0.6981. 



When /i = 1/2 the relaxation (nndging) parameter is too large for the 
approximating solntion to converge to the reference solntion, and when /z = 
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1/5 it is too small. However, the intermediate value /i = 1/3 works with the 
error represented by ||M(t) — v{t)\\v falling below 10 by t = 13 417.8. Note 
that, since the double-precision floating-point numbers used to represent the 
Fourier modes of u and v on the computer have only 15 digits of precision, 
we don’t expect convergence of \\u(t) — v(t)\\v to exact zero over time. 

For /i = 1/4 the error falls below 10“^° at T = 12 327.1, however, it 
rises again and it is not clear whether after T = 23 463.9 the error hnally 
stays below 10“^° or not. The value /r = 2/9 shows an even more irregular 
pattern where ||n(f) —v{t)\\v exhibits a period of decay followed by a period 
of growth that covers six orders of magnitude. Fortunately, most of our 
parameter choices avoid these borderline cases and the corresponding error 
either converges towards zero and stays below 10“^*^ or shows few signs of 
converging and stays well above 10“^°. 

To determine the values of h and /i for which it is possible to recover the 
reference solution to within numerical roundoff error we £x e > 0 and dehne 

Tmax = sup{f e [0,T] : ||n(f) - u{t)\\y > e} 

and 

Tmin = inf {f e [0,T] : ||n(f) - u{t)\\y < e} 

Let Tjnax = 0, when the supremum is over the empty set, and T^in = oo when 
the inhmum is over the empty set. When ||n(T) — u{T)\\v > e further set 
Tmax = C )0 to eusure Tma x > Tmin. Inspired by Figure [3] we also define 

1 fT 

£avg = 7^ —^ / ||n(t) - u{t)\\vdt 

I 10 Jtq 

and take e = 10“^°, T = 25 000 and Tq = 2T/3 for our numerics. 

Table [T] shows the results of our computational experiments. Runs with 
K = 8 were also performed, however, no value of fi yielded a hnite value for 
Tmax or Tinin or eveu an approximation for which the error \\u{t) — v{t)\\v fell 
below 10“^. We conclude that iF = 9 is the minimal resolution for which 
there exists a fi such that the error tends toward zero. At this minimal 
resolution only a narrow range of values for /i near 1/3, result in an error 
which falls below 10“^°. When K = 10, there is a much greater range of 
corresponding values for /x that work well. In fact, when iF = 10 all values of 
H between 1/6 and 1/2 led to corresponding approximations v{t) such that 
^avg ~ 3 X 10“^^. Since the double-precision floating-point numbers used to 
represent the Fourier modes of u and v on the computer have only 15 digits of 
precision, the fact that the error can approach 10“^^ is remarkable. We again 
note, as is consistent with the analysis in [1], our numerical experiments do 
not recover the reference solution if /x is too small or too large. 
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Table 1: Data assimilation using 1^. 


h 

T 

mm 

K = 9 
T 

max 

'^avg 


T 

mm 

K = 10 

T 

max 

•^avg 

0.0625 

oo 

oo 

1.1 


oo 

oo 

8.8 X 10-1 

0.125 

oo 

oo 

3.9 X 

10-1 

13754.5 

17094.3 

4.1 X 10-12 

0.154 

oo 

oo 

2.4 X 

10-1 

4331.1 

4432.6 

5.4 X 10 - 1 ^ 

0.167 

oo 

oo 

2.4 X 

10-1 

3965.0 

4187.7 

3.1 X 10 - 1 ^ 

0.182 

oo 

oo 

1.5 X 

10-1 

3320.1 

3320.1 

2.8 X 10 - 1 ^ 

0.2 

oo 

oo 

3.5 X 

10-2 

2825.2 

2825.2 

2.8 X 10 - 1 ^ 

0.222 

oo 

oo 

1.9 X 

10-5 

2870.1 

2870.1 

2.5 X 10 - 1 ^ 

0.25 

12327.1 

23463.9 

3.7 X 

10-10 

2701.0 

2701.0 

2.5 X 10 - 1 ^ 

0.286 

12275.6 

13466.2 

1.3 X 

10-12 

2581.4 

2598.6 

2.1 X 10 - 1 ^ 

0.333 

13417.8 

13417.8 

1.4 X 

10-12 

2601.8 

2601.8 

2.0 X 10 - 1 ^ 

0.4 

oo 

oo 

1.3 X 

10-1 

3008.4 

3008.4 

2.1 X 10 - 1 ^ 

0.5 

oo 

oo 

3.2 X 

10-1 

4564.3 

4564.3 

2.4 X 10 - 1 ^ 

0.6 

oo 

oo 

6.7 X 

10-1 

8604.5 

9964.6 

5.7 X 10 - 1 ^ 

0.7 

oo 

oo 

1.5 


oo 

oo 

5.3 X 10-2 


Figure 4: Tmax versus /r for e = 10 T = 25000 and K = 10. 



Next consider the family of smoothed interpolants 1^ for different values of 
rj. Figure IHplots Tma x versus /i. When rj is near 0.7 we hnd that values of /x be¬ 
tween 1/4 and 64 all lead to approximations such that \\u{t) —v(t)\\v < 10“^° 
for large enough T. Thus, smoothing with rj near 0.7 leads to a signihcantly 
wider range of values for /i such that the data assimilation algorithm can be 
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used to recover the reference solution. Note when rj = 2 and K = 10 that 
no values of /i led to the convergence of the approximating solution to the 
reference solution over time. 

Having found good values for rj we continue our numerical study by fixing 
rj = 0.7 and varying /i for different resolutions h = L/K, where K = 8, 9 
and 10. As before e = 10“^°, T = 25 000 and Tq = 2T/3. The computa¬ 
tional results given in Table [2] show that observational measurements with a 
resolution given by A' = 8 can now lead to an approximate solution which 
converges to the reference solution over time. Moreover, the accuracy of the 
resulting approximations also improve compared to the non-smoothed case. 
Note that we have omitted reporting in Table [2] since in all cases T^in 
was equal or nearly equal to T^ax- 


Table 2: Data assimilation using where rj = 0.7. 


A 

K = 8 

T £ 

K = 9 

T £ 

a: = 10 

T £ 

0.25 

oo 5.1 X 10-^ 

CX) 9.9 X 10-^ 

2834.2 2.4 x lO'^^ 

0.5 

oo 6.5 X 10“^ 

2570.1 2.4 X 10-1^ 

1534.2 1.7 X 10-1^ 

1 

2817.8 1.1 X 10-^3 

1686.5 1.8 X 10-1^ 

1112.7 1.6 X 10-1^ 

2 

2527.1 2.6 X 10-^^ 

1232.6 1.7 X 10-1^ 

897.5 1.6 X 10-1^ 

4 

2013.1 2.1 X 10-^^ 

1092.4 1.7 X 10-1^ 

718.9 1.6 X 10-1^ 

8 

2191.6 2.1 X 10-^^ 

1124.8 1.7 X 10-1^ 

717.8 1.6 X 10-1^ 

16 

4137.3 4.7 X 10-^^ 

1360.4 1.7 X 10-1^ 

769.7 1.6 X 10-1^ 

32 

oo 3.0 X 10“^ 

2752.9 2.1 X 10-1^ 

1284.2 1.7 X 10-1^ 

64 

cx) 1.4 

oo 4.9 X 10“^ 

2848.6 2.1 X 10-1^ 

128 

oo 2.6 

oo 1.9 

oo 1.1 


Figure |5] plots the data in Table |2j From this figure it is clear that the 
data assimilation algorithm given by f[T5|l with Ih replaced by Ih works well 
when T] = 0.7 for a wide range of values of the relaxation parameter fi. 


4 Numerical Methods 

All fluid dynamics simulations presented in this paper were performed using 
a new parallel code written for the NVIDIA Compute Unified Device Archi¬ 
tecture in the CUDA C programming language [TB] which was developed on 
desktops at the University of Nevada Reno and run on the Big Red II Cray 
XE6/XK7 supercomputer at Indiana University. Computations were made 
in the stream function formulation using a fully-dealiased spectral Galerkin 
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Figure 5: T^ax versus /i for e = 10 T = 25000 and r] = 0.7. 



method. Time steps were performed using a split-Euler method in which the 
linear term was integrated exactly and the non-linear terms were integrated 
using forward differences. 

Specifically we set AT = curl u and compute the reference solution using 
the stream function formulation 

^^-z/A 2 t + /3(T) =curl/, (16) 

where 

/3(T) = J(T, AT) = T,AT, - T^^AT, 

= ((^x) ~ (^J/) )xy ~ {^x'^y)xx + x^y)yy 

Similarly, set AT = curl v and compute the approximating solution using 

^ - ^.A^T + /3(T) = curl / - /r(i?,(T) - /?,(T)), (17) 

where i?/i(T) = curl Po-J^(curl~^ AT). Note that Rh- —>■ V~^. 

Following the 2/3 dealiasing rule applied to 512x512 sized discrete Fourier 
transforms we set JC = {—341,..., 341}^ and approximate 

T ^ Tfc T ^ 5] Tfc P"'" and curl / = ^ 

k£lC k&K. k£K 

Substituting these approximations into flT6|l and flTTl) and projecting onto the 
Fourier modes with wave numbers k & K, yields the Galerkin truncations 

- T^^kk'^ + /3(4')fc = h 
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and 


= g,- 

The corresponding nnmerical scheme for the reference solntion is 
+ At) ^ 

and for the approximating solntion is 

Mt + At) ^ + ^{kHt))k + fiRhi^ - ^)fc)} 

- 47^fc(l 

At the discrete level it is still the case that only nodal-point observational 
measnrements of the reference solntion are nsed to constrnct the approximat¬ 
ing solntion. Moreover, since both solntions are integrated nsing the same 
nnmerical methods, we may think of T as an nnknown reference solntion 
that evolves according to a known discrete dynamical system, and the soln¬ 
tion represented by $ as an approximation generated by data assimilation 
according to the exact same discrete dynamics. Thns, even thongh onr nn¬ 
merical schemes are only first order in time, we consider onr nnmerical ex¬ 
periments to simnlate data assimilation in the absence of both measnrement 
and model errors. 

We take the time step to be At = 1/2048 which is enongh to ensnre that 
the CFL condition 

NAt ,, , ,, , 

snp {|Mi(a;)| \u 2 {x)\} < 0.1608 -C 1 

2T a:gS7 

is satished for the reference solntion over the entire rnn. Note that as fi 
gets larger the data assimilation eqnations flT5|l become stiffer. Therefore, 
the time step At was also chosen small enongh to ensnre the stability of the 
conpled nnmerical scheme for compnting the approximating solntion. 

Onr nnmerical software has been optimized so that it rnns entirely on the 
CUBA hardware with zero memory copies and fonr fast Fonrier transforms 
per time step. Note that fonr transforms is the minimal nnmber for the two- 
dimensional Navier-Stokes eqnations, see Basdevant [2] for farther remarks 
and analogons optimizations when compnting the three-dimensional Navier- 
Stokes eqnations. As device memory on the CUBA hardware is relatively 
scarce we also minimized onr storage reqnirements. Storage reqnirements 
consist of fonr n x n donble-precision scalar arrays: one for T, another for 
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Figure 6 : Data Flow for Computing the Non-linear Term 
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$ and two temporary arrays Ti and T 2 . Figure | 6 ] shows the data flow when 
computing the non-linear term. The first line represent the contents of Ti, 
the second represents T 2 and the arrows represent computational kernels. 

When using 512 x 512 FFTs, our code achieves approximately 1062 time 
steps per second. In particular, the computational speed running on an 
NVIDIA Tesla K20 GPU was found to be roughly 37 times faster than equiv¬ 
alent CPU code running on single AMD Opteron 6212 core and 11.5 times 
faster when compared to running on 32 CPU cores. Correctness of operation 
was verihed using the Navier-Stokes solver described in [El and [IB]. 


5 Conclusion 

As is consistent with the analytical bound ([ 8 ]) and related discussion in [ 1 ], the 
numerical results given in Table [T] and Figure [5] show that the approximating 
solution does not converge to the reference solution when /r is either too small 
or too big. At the same time, provided the resolution h is fine enough and 
77 ~ 0.7, there is a wide range of good values for /i when using the smoothed 
interpolant observable Ih- In particular, when h = L/10 ~ 0.6283, the data 
assimilation algorithm ([7]) performs similarly for values of /r between 0.5 and 
32. Note, however, that smaller values of /i are computationally preferable 
because of stiffness considerations. 

We remark that our numerical experiments have been conducted using 
exact error-free measurements and exact model dynamics and that in the 
presence of measurement and model errors we don’t expect a similarly wide 
range of good values for p. In fact, preliminary computations show when is 
noise added to the system that there exists a unique optimal value for fi re¬ 
flecting the tradeoff between measurement and model errors. Theoretically, if 
the dynamics represented by J^{u) in ([T]) are linear, then /r can be seen as the 
parameter of a linear Kalman filter, see for example Majda and Harlim [15] . 
and there exists an analytically derived optimal value for /i which represents 
the tradeoff between measurement and model errors. In the fully non-linear 
case studied here, the fact that there is a wide range of good constant values 
for fi in the absence of measurement and model errors suggests, provided the 
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resolution h is fine enough, that can be further optimized as if the model 
was linear. 

We now compare the coarsest resolution h = L/K that works for the 
data assimilation experiments presented here when rj = 0.7 to the number 
of numerically determining modes found in [18]. Under the same physical 
parameters and forcing, the minimum number of Fourier modes needed by 
([2D was 

Uc = card(P5) = 80 

where 

VR = {e^'^-^-QKkl + kl^R^}. 

In this paper we show the minimum of nodal measurements needed are 

N = K‘^ = U 

which by the Nyquist-Shannon sampling theorem may be represented by the 
Fourier modes 


■^K = { : 0 < max(|fci|, \k 2 \) < Kj2 }. 

To compare these two results we note that T>r represents a circle in Fourier 
space while Mk represents a square. Let i?min = 5 and iFmin = 8. If the 
resolution requirements of algorithm ([2D are comparable to ([7D, we would 
expect that Mk ^ , ^ ^ would imply R > i?min and that C Mr would 

imply K > iFmin- This is supported by our results. If TVg C Vr then 

R>^^ 5.65 > 5 = R^,^. 

Similarly, if C Mr then 

K>2-5 = 10>8 = 

Thus, even though our nodal observations possess the problems of aliasing 
and high-frequency spill over, these problems can be mediated with appro¬ 
priate smoothing. The resulting resolution K needed for the approximating 
solution to converge to the exact solution is then about the same as suggested 
by the number of numerically determining modes. 
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